function [best] = linear_svm(data)
	%addpath(genpath('../liblinear-1.92')); <- move into setup.m
	
	% Recommended by libsvm
	Cs = [2^-5, 2^-3, 2^-1, 2^1, 2^3, 2^5, 2^7, 2^9, 2^11, 2^13, 2^15];
	C = length(Cs);

	results = {};
	results.TR.acc = zeros(C,1);
	results.TR.p    = cell(C,1);
	results.TR.dv   = cell(C,1);
	results.V.acc  = zeros(C,1);
	results.V.p     = cell(C,1);
	results.V.dv    = cell(C,1);
	results.TE.acc = zeros(C,1);
	results.TE.p    = cell(C,1);
	results.TE.dv   = cell(C,1);
	results.models  = cell(C,1);
	results.libsvmmodel = cell(C,1);
	
	% class weights
	pos  = sum(data.YTR == 1);
	neg = sum(data.YTR == -1);
	
	% [data.XTR, data.YTR] = rebalance(data.XTR, data.YTR);
	% [data.XTV, data.YTV] = rebalance(data.XTV, data.YTV);
	% [data.XTE, data.YTE] = rebalance(data.XTE, data.YTE);

	for c = 1:C
		%['-c ' num2str(Cs(c)) ' -B 1 -q']);%
		model = train(data.YTR, sparse(data.XTR), ['-c ' num2str(Cs(c)) ' -B 1 -w1 ' num2str(neg/(neg+pos)) ' -w-1 ' num2str(pos/(neg+pos)) ' -q']);
		[ptr, ~, dtr] = predict(data.YTR, sparse(data.XTR), model); 
		[ptv, ~, dtv] = predict(data.YTV, sparse(data.XTV), model); 
		[pte, ~, dte] = predict(data.YTE, sparse(data.XTE), model); 

		results.TR.acc(c) = sum(ptr == data.YTR)/length(data.YTR);
		results.TR.p{c}   = ptr;
		results.TR.dv{c}  = dtr;
		results.V.acc(c)  = sum(ptv == data.YTV)/length(data.YTV);
		results.V.p{c}    = ptv;
		results.V.dv{c}   = dtv;
		results.TE.acc(c) = sum(pte == data.YTE)/length(data.YTE);
		results.TE.p{c}   = pte;
		results.TE.dv{c}  = dte;
		results.models{c} = model.w';
		results.libsvmmodel{c} = model;
		
	end
	
	[M,CIX]  = max(results.V.acc);
	
	% best = {};
	% best.TR.acc = results.TR.acc(CIX);
	% best.TR.p   = results.TR.p{CIX};
	% best.TR.dv  = results.TR.dv{CIX};
	% best.V.acc  = results.V.acc(CIX);
	% best.V.p    = results.V.p{CIX};
	% best.V.dv   = results.V.dv{CIX};
	% best.TE.acc = results.TE.acc(CIX);
	% best.TE.p   = results.TE.p{CIX};
	% best.TE.dv  = results.TE.dv{CIX};
	% best.model  = results.models{CIX};
	% best.libsvmmodel = results.libsvmmodel{CIX};
	
	% % train on training and validation
	% ['-c ' num2str(Cs(CIX)) ' -B 1 -q']);
	model = train([data.YTR; data.YTV], sparse([data.XTR; data.XTV]), ['-c ' num2str(Cs(CIX)) ' -B 1 -w1 ' num2str(neg/(neg+pos)) ' -w-1 ' num2str(pos/(neg+pos)) ' -q']);
	%['-c 1 -B 1 -w1 ' num2str(neg/(neg+pos)) ' -w-1 ' num2str(pos/(neg+pos)) ' -q']);
	[ptr, ~, dtr] = predict(data.YTR, sparse(data.XTR), model); 
	[ptv, ~, dtv] = predict(data.YTV, sparse(data.XTV), model); 
	[pte, ~, dte] = predict(data.YTE, sparse(data.XTE), model); 
	
	best = {};
	best.TR.acc = sum(ptr == data.YTR)/length(data.YTR);
	best.TR.p   = ptr;
	best.TR.dv  = dtr;
	best.V.acc  = sum(ptv == data.YTV)/length(data.YTV);
	best.V.p    = ptv;
	best.V.dv   = dtv;
	best.TE.acc = sum(pte == data.YTE)/length(data.YTE);
	best.TE.p   = pte;
	best.TE.dv  = dte;
	best.model  = model.w';
	best.libsvmmodel = model;
	best.C      = Cs(CIX);
	
end
	
	